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Anomalous diffusion of pions at RHIC 
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After pointing out the difference between normal and anomalous diffusion, we consider a hadron 
resonance cascade (HRC) model simulation for particle emission at RHIC and point out that rescat- 
tering in an expanding hadron resonance gas leads to a heavy tail in the source distribution. The 
results are compared to recent PHENIX measurements of the tail of the particle emitting source in 
Au+Au collisions at RHIC. In this context, we show how can one distinguish experimentally the 
anomalous diffusion of hadrons from a second order QCD phase transition. 
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" It does not make any difference 
how beautiful your guess is. 
It does not make any difference 
how smart you are, 

who made the guess, or what his name is - 
if it disagrees with experiment, it is wrong. " 

/R.P. Feynman/ 



I. INTRODUCTION 

Various new techniques are being developed in a vi- 
brant and inspiring, sometimes challenging and puzzling 
sub-field, named recently by Lednicky [lj as correlation 
femtoscopy. A series of inspiring and stimulating recent 
reviews @, S 0, IE @ re-connected femtoscopy to the 
search of new phases of QCD as well as other new, some- 
times unexpected, sometimes puzzling expectations and 
observations. One of these new experimental results has 
been a recent PHENIX measurement of source images in 
Au+Au collisions at ^Jsnn = 200 GeV that found evi- 
dence for a non-Gaussian structure and a heavy (larger 
than Gaussian) tail in the source distribution for pions 
Q. This points to an interesting new direction, going 
beyond investigating only the means and the variances 
of the source distributions on the scales of femtometers. 

Resonance decays are known to be able to produce 
long, non-Gaussian tails, because some of the resonances 
have large decay times as compared to the characteris- 
tic 4-5 fm/c source sizes extracted from interferomctric 
measurements in Au+Au collisions at RHIC. In fact, a 
smaller and smaller fraction of resonances have larger 
and larger life-times, and this effect was considered first 
by Bialas Q who argued that even a cut power-law cor- 
relation functions may appear, due to the fact that the 
decay times of pion producing resonances have a broad 
probability distribution. Another possible conventional 



source for such a heavy tail of the source of pions might 
be the elastic rescattering of the produced hadrons: as 
the hadron gas expands, the system becomes cooler and 
more and more diluted, hence the mean free path be- 
comes larger and larger. In many hydrodynamical calcu- 
lations, an idealized freeze-out process is assumed, when 
the mean free path suddenly jumps from (the hydro- 
dynamic limit) to infinity (the freeze-out limit). More 
realistically, the mean free path diverges to infinity in 
a finite time interval, and rescattering in a time depen- 
dent mean free path system is known to lead to new 
phenomena, which has been studied in great detail un- 
der the name of anomalous diffusion in other branches 
of physics. One of the experimentally observed charac- 
teristics of such an anomalous diffusion pattern was the 
appearance of approximately power-law shaped tails in 
the coordinate space distributions, which is to be con- 
trasted to the Gaussian, strongly decaying tails observed 
in normal diffusion or in Brownian motion. We shall ex- 
plore the phenomena of anomalous diffusion in the con- 
text of heavy ion physics, first pointing out its general 
mathematical properties and its relations to Levy source 
distributions, following the review of Klafter and Metzler 
on anomalous diffusion Q . Then we shall explore anoma- 
lous diffusion of pions, kaons and protons in Au+Au col- 
lisions at y/SNN — 200 GeV colliding energies using the 
simplest possible tool, namely the conventional Hadronic 
Resonance Cascade (HRC) model of Tom Humanic [10j. 
First we compare the results of this HRC simulation with 
PHENIX data. Then we investigate the sensitivity of 
the characteristics of the simulation for various experi- 
mentally available controls, like the selection of the cen- 
trality class of the events or the selection of a transverse 
momentum region of the pair. Finally we summarize and 
conclude. As a starting point, let us consider a simpli- 
fied version of the presentation in ref. [9(, to review the 
equations of normal and anomalous diffusions in the same 
framework, based on a master equation approach. 
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II. NORMAL DIFFUSION 

Normal diffusion corresponds to a physical process 
when one investigates the motion of a test particle in 
a medium, which has some grid like structure with fixed 
(say, Ax) lattice constant, (equidistant cells) and jumps 
(of a fixed frequency) are allowed to nearest neighbor- 
ing cells. Let us for simplicity consider a one dimen- 
sional, normal diffusion process in the framework of mas- 
ter equation approach. The time is denoted by t, A test 
particle is initially located at a given cell, denoted by 
j = 0. Different cells are indexed by the integer j £ Z. 
The probability distribution of finding the particle in cell 
j at time t is denoted by Wj (t) . (Note that Wj (t) can be 
considered as analogous to the particle emission function 
S(x,t), that we frequently encounter in particle inter- 
ferometry and femtoscopy.) Suppose that particles may 
jump from the given cell to its nearest neighbors, ran- 
domly up and down, at certain regular but small time in- 
tervals At. The continuum limit means At — > 0, Ax — > 
in such a way that the mean free path (and also the mean 
collision time) remains constant. 

In this situation, the following master equation drives 
the time evolution in this material: 



1, 



1, 



Wj(t + At) = -Wj-i® + ^W J+1 (t). 



(1) 



We rewrite this discretized form to a continuous form 
by introducing the continuous coordinate variable x. We 
assume that the Ax step size is infinitesimally small com- 
pared to the overall length-scale of the medium, and that 
the time interval between the subsequent jumps is small 
as compared to the time duration corresponding to the 
observation of the diffusion process. These assumptions 
yield the following leading order Taylor expansions: 

W J (t + At) = W 1 (t) + At^+0(At 2 ), (2) 

W J±l{ ^W { x,t)±Ax d -l + ^f d -^ + 0(Ax% 

(3) 

and these expressions can be combined to derive the con- 
tinuum form of the diffusion equation 



dW d 2 , , 

- W =K 1 —W(x,t). 



(4) 



All properties of this matter are characterized by the dif- 

A 2 

fusion constant K\ — limAa;-+o,At->o ^St , which corre- 
sponds to the mean squared displacement per unit time. 

The above (normal) diffusion equation can be solved 
easily by introducing the Fourier-transform 



W(k,t) 



dx exp(ikx)W(x, t), 



(5) 



that leads to the momentum-space diffusion equation 

^ = - Kl k 2 W(k,t). (6) 



The solution of this equation is a Gaussian function of 
k that can be converted back to the coordinate-space 
distribution to yield the solution of the normal diffusion 
equation: 



W(x,t) 



1 



exp 



(7) 



corresponding to normal or Gaussian diffusion of test par- 
ticles. The initial condition corresponding to this solu- 
tion is indeed 



W (x) = \\mW{x,t) = S(x) , 



(8) 



corresponding to a localized package inserted in a ma- 
terial with homogeneous, time independent properties. 
Hence in a homogeneous and time independent medium, 
the diffusion of point-like initial source yields a Gaus- 
sian source density distribution, and the mean square of 
this Gaussian, R 2 = 1K\t increases linearly with increas- 
ing time. This behavior will be contrasted to the more 
general case, when both the jump lengths and the jump 
frequencies have a continuous probability distribution. 



III. ANOMALOUS DIFFUSION 

We again follow Metzler and Klafter in deriving 
anomalous diffusion from the so-called continuous time 
random walk models. Suppose that the jump length and 
jump time have the probability distribution ip(x, t). Then 
the jump length distribution X(x) and the waiting time 
or jump time distribution w(t) reads as 



X(x) 
w(t) 



4>(x, t) dt, 
ip(x, t) Ax. 



(9) 
(10) 



Suppose that an external force F(x) also influences the 
motion of the test particles. These considerations lead to 
the generalized or anomalous diffusion equation, which 
is a kind of a generalized Fokker-Planck equation for the 
phase-space distribution, W(x,v,t): 



dW 
~dt 



dW F(x) dW 



dx 



dv 



l-a' 



W, (11) 



which contains fractional derivatives and other subtleties 
that go well beyond the scope of this presentation. For 
the detailed definition and explanation of this fractional 
Fokker-Planck equation we refer to Q. What is impor- 
tant, however, that if the waiting time distribution has 
a Poissonian shape, the exact solution of this fractional 
Fokker-Planck equation has a simple form in the momen- 
tum space representation: 



W(k,t) = eM~tK a \k\ a ). 



(12) 
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This form is the well known characteristic function 
(Fourier-transform) of Levy stable source distribu- 
tions [HI, [13, EH E3 • Here the parameter a stands for the 
Levy index of stability, in general < a < 2 for Levy sta- 
ble source distributions, and parameter K is an anoma- 

A 2 

lous diffusion constant, K — limAz^o.Ai^o 2At ^ /Q ■ 

Note that Levy stable distributions were introduced to 
particle interferometry studies recently in ref . [l5j , based 
on general, mathematical arguments like generalized cen- 
tral limit theorems. Anomalous diffusion is a specific ex- 
ample of a physical process that under certain conditions 
detailed in [9j satisfies such generalized central limit theo- 
rems: one more rescattering in the diffusion process does 
not change the limiting behavior of the source distribu- 
tion. 



IV. LEVY STABLE LAWS AND ANOMALOUS 
DIFFUSION 

The edification from the previous section is that rescat- 
tering in a system with a time dependent mean free path 
under certain conditions (e.g. Poissonian waiting time 
distributions) leads to a random Levy walk (or anoma- 
lous diffusion), instead of Brownian motion (or normal, 
Gaussian diffusion). In case of anomalous diffusion or 
Levy walks, the scale parameter grows with time as 
R a tx t, in contrast to normal diffusion, that corresponds 
to the a = 2, R 2 oc t special case. 

The difference between normal and anomalous diffu- 
sion is illustrated in Fig.[TJ reproduced from ref. @. 

Although the mean displacement or the variance of 
these distributions might diverge, one additional step in 
the anomalous diffusion does not change the limiting be- 
havior of the process. 

It is also interesting to note that the tail of anomalous 
diffusions has a power-law structure in coordinate space, 
S(r) oc for r 3> R, where d is the number of 

spatial dimensions (d = 3 in our world), and a is the same 
Levy index of stability as before. This asymptotics is 
yet another important property of the Levy stable source 
distributions. 

Nature often violates Gaussian universality, mirrored 
in experimental results which do not follow Gaussian pre- 
dictions Q. The evidence for a non-Gaussian behavior in 
multivariate and univariate source distributions has been 
observed recently also in high energy heavy ion collisions. 
In fact, the first observation of a non-Gaussian corre- 
lation function in Au+Au collisions at RHIC has been 
made by the STAR collaboration that utilized an Edge- 
worth expansion (Lsj and quantified the deviation from 
the Gaussian structure of the particle emitting source in 
terms of non-vanishing fourth order cumulant moments 
of the squared Fourier-transformed source distribution. 
More recently, the PHENIX Collaboration has applied 
the imaging method of Brown and Danielewicz [171 ] to 
reconstruct the two-particle relative coordinate distribu- 
tion Q. In this paper, PHENIX observed a clear de- 




FIG. 1: Simple illustration of the qualitative properties of 
normal diffusion (bottom) and anomalous diffusion (top). In 
the latter case, large jumps separate various more local parts 
of the trajectory. The size of the covered region is larger in 
case of the anomalous diffusion than in case of the Gaussian, 
normal diffusion. From ref. JQj. 



viation from a Gaussian structure, and pointed to the 
appearance of a heavy tail in the two-particle relative 
coordinate distribution. In the subsequent part of this 
manuscript, we investigate if these PHENIX source func- 
tion data can be reproduced with the help of Monte-Carlo 
models that incorporate the concept of anomalous diffu- 
sion and that are able to describe the more standard 
observable like single particle spectra and the three di- 
mensional Gaussian fit parameters, i? S ido, -Rout and i?i on g 
of the measured two-pion correlation functions in high 
energy heavy ion collisions. 



V. MONTE-CARLO SIMULATIONS OF 
ANOMALOUS DIFFUSION OF PIONS 

Heavy ion collisions produce thousands of particles in 
a single high energy nucleus- nucleus collision. The bulk 
of the particle production, i.e. the momentum distribu- 
tion and the correlation patterns of 99 % of the particles 
are best described by hydrodynamical, or hydrodynam- 
ically inspired models, see ref. 0] for a recent review. 
Now we are interested in the tails of particle production, 
which might indicate a deviation from the hydrodynam- 
ical behavior. Hence our attention is turned to Monte- 
Carlo (MC) simulations. Two Monte-Carlo models, the 
Hadronic (or Humanic) Resonance Cascade Model [l8| . 
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(HRC) and the AMPT model of Zhang^Ko, Li and 
Lin [l9j] have also demonstrated [13, H3, 21\ their abil- 



ity to describe single particle spectra, elliptic flow, and 
HBT correlation measurements in Au+Au collisions with 
^/JW^ = 130 and 200 GeV at RHIC. 



A. Selection criteria - comparison with data 

The selection of the MC model was based on the follow- 
ing criteria. We looked for a conventional hadronic cas- 
cade model that describes single particle spectra, elliptic 
flow data, and HBT data (without any puzzles), hence 
yields a good description of the hadronic final state, i.e. 
it yields an acceptable model of S(r). 

It also has to be well documented and easy to use, has 
to work at CERN SPS as well as at RHIC energies, con- 
tain the most important short and long lived resonances 
e.g. uj, r\ and r( (hence can be used as a realistic model 
for halo that appears from the decay products of these 
long-lived resonances). It is important to require that 
the simulation takes into account the rescattering among 
the hadrons due to possible development of a power-law 
tail from the anomalous diffusion, discussed in the earliei 
sections. 

Finally we utilized the Hadronic Rescattering 
Model [18], as it satisfied all the criteria listed above. 
We expect that similar results are obtained with the 
AMPT model, but we did not yet investigated the 
predictions of this code in detail. The AMPT model 
is a multi-phase transport model, while the HRC is a 
simpler hadronic resonance cascade. When looking for 
new effects related to the anomalous diffusion in a time 
dependent mean free path environment, we have opted 
for the simplest possible choice - so that the results 
then can be uniquely related to the hadronic final state 
effects. 



B. Self-consistency criteria - subdivision test 
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FIG. 2: Pion pt distributions for I = 1 and 1 = 5, from ref. [Io| ]. 
The insensitivity of the transverse momentum distribution to 
the number of subdivisions indicates that the scaling proper- 
ties of the transport equations are satisfied by the HRC code. 
Similar subdivision tests are satisfied by the HRC code for 
the elliptic flow as a function of the transverse momentum 
v2(p t ), the elliptic flow as a function of the longitudinal an- 
gular variable r\ as V2(r/), see ref. [13] for further details. 



GeV. Thus this HRC code passed not only the test of 
reproducing the important global observables, but also 
the subdivision test, which suggests that there are no 
significant programming artifacts, or causality violating 
terms in this Monte Carlo simulation code. Fig. [5] (taken 
from ref. [13]) shows a comparisons of the I = 1 and 
the / = 5 cases utilizing the HRC model to describe the 
charge averaged pion spectra. Similar plots were pub- 
lished for the HBT radius parameters and the transverse 
momentum and pseudorapidity rj dependence of the el- 
liptic flow. See ref. [13] for further details and for the 
comparison plots with RHIC Au+Au data. 



The self-consistency check of subdivision invariance is 
based on the invariance of the Boltzmann equation, the 
basis of the Monte Carlo particle-scattering calculations, 
for a simultaneous decrease of the scattering cross sec- 
tions by some factor Z, and an increase of the particle 
density by the same factor of Z[22j|. As I becomes suf- 
ficiently large, non-causal artifacts become insignificant. 
The HRC rescattering calculations have been tested by 
comparing pion observables for no-subdivision, i.e. 1=1, 
with subdivision of I = 5. Descriptions of how the ob- 
servables are extracted from the rescattering calculation 
are given elsewhere [20| . 

All calculations were carried out for an impact parame- 
ter of 8 fm, to simulate semi-central collisions, that result 
in significant elliptic flow. Such HRC results [13] repro- 
duced reasonably well the RHIC data on the correspond- 
ing observables in Au+Au collisions at ^Jsnn = 200 



C. Self-consistency criteria - comparison with 
exact hydrodynamical results 

Motivated by the success of the Buda-Lund parame- 
terization [23|, [24| of the hadronic final state, a search 
started to find parametric, but time dependent, exact so- 
lutions of hydrodynamics, that lead to Buda-Lund type 
of nearly Gaussian freeze-out distributions. Surprisingly 
large classes of exact, parametric solutions of hydrody- 
namics have been recently found this way, both in the 
non-relativistic [2a . [2 61 [27| as well as in the relativistic 
[H [H [13, HI [llTsJkinematic domain. 

These generalized, Buda-Lund type parametric, hydro- 
dynamical calculations have a built-in self-quenching ef- 
fect, as analyzed in detail in ref. [27| • For example, all 
the HBT radii stop to evolve in time, they approach a 
direction independent constant, and expansion velocities 
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FIG. 3: Time evolution of an expanding ellipsoid with Gaus- 
sian density profile and homogeneous temperature profile in 
Buda-Lund type of exact solutions of non-relativistic hydro- 
dynamics. (X, Y, Z) stand for the principal axis of this ex- 
panding ellipsoid, after an initial acceleration, these scales 
evolve linearly with time (top). The corresponding HBT radii, 
(Rx,R y ,Rz) approach a direction independent constant (be- 
low top). The expansion velocities of the principal axis of 
the exploding ellipsoid, (V x ,V y ,V z ) tend to direction depen- 
dent constants (above bottom panel) and similarly, the slope 
parameters of the single particle spectra in the principal di- 
rections, (T x ,T y ,T z ) also tend to direction independent con- 
stants (bottom panel). Based on ref. [27t ] . 
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FIG. 4: Time evolution of V2, HBT, and mr slope parameters 
from HRC for RHIC Au+Au, from ref tlC|. The HRC code 
shows the same self-quenching properties as the Buda-Lund 
type of exact solutions of non-relativistic hydrodynamics. 



tend to direction dependent constants, although the sys- 
tem keeps on expanding (via rescattering process or hy- 
drodynamical evolution). Elliptic flow also freezes out at 
the same time when the spectra (slopes) stop to evolve in 
time. These general, qualitative properties of exact para- 
metric hydrodynamic solutions are illustrated in Fig. [3] 
In ref. [27| we have shown that this self-quenching is a 
general property of a large class of exact analytic hydro- 
dynamic solutions, which is independent of the particular 
initial conditions - a beautiful exact result. This property 
of the exact non-relativistic ellipsoidal hydrodynamic so- 
lutions is beautifully reproduced in the Hadronic Reso- 
nance Model, as shown on Fig. 2J 
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D. More than hydro - tails of particle production 

As summarized above the HRC resonance cascade 
model describes well the observables in heavy ion col- 
lisions at both CERN SPS and at RHIC energies, and it 
passes the subdivision test, and yields such a time evo- 
lution of the observables, which is similar to the analyti- 
cally obtained, exact hydrodynamical asymptotic behav- 
ior of these observables. So HRC looks to be a reliable 
model for the production of the bulk of the particles. In 
the subsequent part, we investigate if the HRC model is 
able to describe also the tails of the particle production 
in y / Siv77 = 200 GeV Au+Au collisions. 

It turns out that this conventional hadronic cascade 
model HRC has a built-in adaptive bin size in the time 
direction, hence there is no built-in cutoff time scale in 
the code. HRC also contains the cascading of the most 
abundant hadrons: p, A, K*, id, n, n' , 4> and A, but it ne- 
glects electrical charge. Therefore we see that this HRC 
model implements rescattering in a time dependent mean 
free path system, and in the previous section we have 
shown that under certain conditions this corresponds to 
a random Levy walk and is signaled by power-law tails 
in the source distribution. Let's see whether the HRC 
simulation results confirm such an expectation or not. 



VI. DETAILED HRC SIMULATION RESULTS 

We generated 48 events with b=4.45 fm and 5730 
events with 12.5 fm corresponding to the 0-20% and 50- 
90% centrality classes of PHENIX events. The gener- 
ated events have average multiplicities of 3400 and 38 
particles, respectively. The mean impact parameter val- 
ues suiting the two centralities were determined from a 
Glauber calculation as in ref. (3~il |. 

We made cuts on the data sample similar to the ones 
in the PHENIX imaging paper Q: 

• 0-20% and 0.2 GeV/c < p t < 0.36 GeV/c, 

• 0-20% and 0.48 GeV/c < p t < 0.6 GeV/c, 

• 40-90% and 0.2 GeV/c < p t < 0.4 GeV/c, 

in addition to the PHENIX geometry cut of —0.5 < y < 
0.5. 

The resulting plots are shown in Figs. |5TI2"2"1 The com- 
plete source function is decomposed to different compo- 
nents, i.e. the S{r) distributions created from pairs with 
both pions being 

• primordial or decay products of resonances that 
have a lifetime less than 20 fm/c; this means here 
direct pions p, A, K* decay products, referred to 
as core, or c pions; 

• decay products of an u> 



• decay products of resonances that have a lifetime 
greater than 25 fm/c; in HRC: A, <!>, n, n' , referred 
to as halo or h pions. 

One of the pions in a pair is from one above category, and 
the other may be from another one: so in this study we 
distinguish (c, c), (c, u>), {c,h), {id, id), {id,h) and {h, h) 
type of possible combinations. 



A. HRC simulations and PHENIX data 

In this subsection we compare the HRC model calcu- 
lations with experimental data from [7j. The kinematic 
cuts correspond to Figs. la,b and Fig. 2 of ref. 

The PHENIX imaged S{r) source distribution coin- 
cides with the HRC simulation result, and both are de- 
scribable with a power-law tail, in the kinematic region 
of 0-20% centrality and 0.2 GeV/c < p t < 0.36 GeV/c, 
for pion pairs, as indicated in Fig.O On this log- log plot, 
the tail of S(r) is approximately linear, so it can well be 
called a heavy tail. Clearly a further refinement of the 
experimental resolution could reveal more details on the 
structure of this tail behavior, but a power-law approx- 
imation is not inconsistent with the currently available 
data. 

In HRC, this power-law tail is due to rescattering 
because of the adaptive time-scale, corresponding to 
anomalous diffusion. This results in Levy distributions, 
in contrast to Gaussian distributions, which clearly fail 
to reproduce the tails of the particle production. The 
Levy index of stability, a, can be determined approxi- 
mately from fits to the tails of these S(r) source func- 
tions: as for Levy sources S{r) oc i n d spatial 
dimensions. In the HRC simulations, as well as in case 
of S(r) reconstructed by PHENIX, d = 3. From Fig. [5J 
the HRC simulations and PHENIX data both indicate 
a Levy index of stability of a « 1.15 ± 0.1, which is 
consistent with the direct Levy fits to PHENIX prelimi- 
nary Coulomb corrected correlation functions in a similar 
kinematic region, as presented in ref. [35} . 

This is a great success of the HRC model, as it seems 
that it has implemented the key effect, the anomalous 
diffusion, reasonably well to the simulation. 

Let us investigate in greater detail if HRC can describe 
more subtle features of the PHENIX data. In Fig. [6j we 
show a comparison with data in a similarly soft domain 
but in the more peripheral centrality class. 

Within errors, the HRC simulations again reproduce 
the PHENIX measured S{r) source functions, even if the 
experimental statistics does not allow for a comparison 
in the interesting region of r > 20 fm. The HRC simula- 
tions, which have better statistics than the experimental 
data, indicate a power-law tail with similar exponent as 
in the case of more central collisions. From the com- 
parison with the previous figure it follows that the HRC 
model describes the centrality dependence of the heavy 
tail in this soft k t region in an acceptable manner. This 




FIG. 5: Pions for 0-20% centrality and 0.2 GeV/c < p t < 
0.36 GeV/c, with the various components of the source. On 
the length scales less than 50 fm, core-core pairs are the most 
abundant ones in the present situation. The next to largest 
contribution is from the (c, uS) pairs, however, their contribu- 
tion of is very small as compared to the core-core pairs in the 
experimentally resolvable region, all the other contributions 
are negligible. A comparison of the full HRC simulation re- 
sult (open diamonds) with PHENIX data 0] (filled circles, 
S(r) exp ) indicates that HRC model simulations are in a good 
agreement with data in this kinematic range. Solid black line 
shows the best Gaussian fit to S(r) in the fm < r < 15 fm 
region, which misses the tail region clearly. 



is yet another feature of the HRC model which deserves 
appreciation. 

Let us now investigate how well can a HRC simula- 
tion describe the transverse momentum dependence of 
the relative coordinate distributions of pion pairs. Such 
a comparison is shown in Fig. [7j In the 0-20% centrality 
and 0.48 GeV/c < p t < 0.60 GeV/c kinematic selec- 
tion, the HRC model simulations are in a disagreement 
with PHENIX data. So the HRC model does not de- 
scribe the observed transverse momentum dependence 
of the PHENIX imaged source distribution in case of 
nearly central collisions. Although the achievements of 
this model are really impressive, in this kinematic region 
the model could be improved or fine-tuned. In fact it 
seems that rescattering effects are too large as compared 
to data with increasing transverse momentum. 

In the subsequent parts, we explore the HRC model 
predictions in greater detail. We focus our attention to 
the tails of particle production, and try to identify within 
the limitations of this model calculation what kind of ex- 
perimental control is available for changing the exponent 
(or, on a log-log plot, the slope parameter) of the tails of 
particle emission in the HRC simulations. This is moti- 
vated by the recent predictions in ref. [36| that suggested 
the existence of a power-law tail in the coordinate space 
distribution at the critical end point of the line of first 
order phase transitions in QCD. At this critical point the 
phase transition is of second order, and the Bose-Einstein 
correlation function has a Levy form, with the Levy in- 



FIG. 6: Same as the previous figure, but for a 50-90% cen- 
trality and 0.2 GeV/c < p t < 0.40 GeV/c kinematic selection. 
HRC simulations are shown with the various components 
of the HRC source resolved. A comparison with PHENIX 
data [jj (filled circles, S(r) exp ) indicates that HRC model 
simulations are in an agreement with facts in this kinematic 
range, too, although within errors, systematic deviations be- 
tween the simulation results and PHENIX data can be ob- 
served in the r > 20 fm region. 
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FIG. 7: Same as the previous figures, but for a 0-20% central- 
ity and 0.48 GeV/c < p t < 0.60 GeV/c kinematic selection. 
HRC simulations are shown with the various components of 
the HRC source resolved. The HRC model simulations are in 
a disagreement with PHENIX data in this kinematic range. 



dex of stability coinciding with the correlation exponent 
7] that is one of the critical exponents, and its value is 
universal, depending only on the universality class of the 
second order phase transition. Hence in a second order 
QCD phase transition the Levy index of stability or the 
power-law exponent of the tails of the particle produc- 
tion becomes independent of the momentum range, cen- 
trality selection, and the particle type. However, in case 
of anomalous diffusion of hadrons the rescattering might 
well be sensitive to the centrality that drives multiplic- 
ity and particle densities, the momentum range that also 
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FIG. 8: Source distribution of HRC simulated pion pairs with 
0.2 GeV/c < p t < 0.4 GeV/c for various centrality classes. 



FIG. 10: Source distribution of kaon pairs with 0.2 GeV/c 
< pt < 0.4 GeV/c for various centrality classes. 
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FIG. 9: Source distribution of HRC simulated pion pairs with 
0.5 GeV/c < p t < 1.0 GeV/c for various centrality classes. 



FIG. 11: Source distribution of kaon pairs with 0.5 GeV/c 
< pt < 1.0 GeV/c for various centrality classes. 



influences how many particles take part in the rescatter- 
ing process and also the particle type, as the number of 
rescatterings is expected to depend on the particle cross 
sections. 



B. Centrality dependence of the HRC source 

The centrality dependence of the HRC simulated 
source distributions of pions with various pair transverse 
momenta is shown in Figs. [5KT31 The centrality depen- 
dence of the tails of particle emissions is found to be 
surprisingly small, as on the log-log plots the tails corre- 
sponding to various centrality selections in Figs. [8lfT3l are 
found to be parallel. The centrality selection, however, 
sensitively influences the region of bulk particle produc- 
tion, corresponding to the change of the scale parameters 
or the effective source sizes in the small r region. The 
slope of the tails in this simulation, however, is indepen- 
dent of centrality for pions at low or higher transverse 



momentum, Figs. [8][9l as well as for low or higher mo- 
mentum kaons, Figs. [TUlfTTl and for protons, Figs. [HHfT51 
Although in case of the proton sources statistical lim- 
itations prevent us from firm conclusion, it seems that 
the power-law exponent of the tails of particle emission 
is rather insensitive to the centrality selection, regardless 
of particle type and pair momentum selection, so central- 
ity is not a sensitive control tool to separate heavy tails 
that arise due to anomalous diffusion from heavy tails 
that appear due to the vicinity of a second order QCD 
phase transition. 



C. Transverse momentum dependence 

Transverse momentum dependence of simulated source 
distributions in various centrality classes is shown in 
Figs. [TUfnll For nearly central collisions and pions, the 
increase of the transverse momentum of the pair reduces 
the size of the bulk production region, as was well known 
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FIG. 12: Source distribution of proton pairs with 0.2 GeV/c 
< pt < 0.4 GeV/c for various centrality classes. 



FIG. 15: Source distribution of pion pairs with 40-80% cen- 
trality, for various p t ranges 
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FIG. 13: Source distribution of proton pairs with 0.5 GeV /c 
< Pt < 1.0 GeV/c for various centrality classes. 



FIG. 16: Source distribution of kaon pairs with 0-20% cen- 
trality, for various pt ranges 
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FIG. 14: Source distribution of pion pairs with 0-20% cen- 
trality, for various pt ranges 



FIG. 17: Source distribution of kaon pairs with 40-80% cen- 
trality, for various p t ranges 
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FIG. 18: Source distribution of proton pairs with 0-20% cen- 
trality, for various pt ranges 



FIG. 19: Source distribution of proton pairs with 40-80% cen- 
trality, for various pt ranges 



from earlier Bose-Einstein correlation measurements and 
theoretical explanations, see e.g. refs. [23|, [24[. Being 
aware of the strong transverse momentum dependences 
of the scale parameter R, which is well shown by the HRC 
simulation, it is rather surprising, that the shape param- 
eter of the tail, the Levy index of stability a is remark- 
ably insensitive to the selected transverse momentum re- 
gions, which follows from the fact that on the log-log plot 
the simulated S(r) tails are parallel to one another, see 
Fig. [2J The same effect is shown in Fig. [TS1 in case 
of peripheral collisions. For kaons, the evolution of the 
source function is less apparent in the same transverse 
momentum regions, this is due to the fact that the effec- 
tive radius parameters, the scales in analytic calculations 
like of refs. [HI, |24j, are found to be predominantly de- 
pending on the transverse mass, m t = \J m 2 + pf of the 
particles. The same variation in the transverse momen- 
tum pt leads to a larger variation in mt for pions, as 
compared to that of kaons, due to the larger value of the 
kaon mass: w 140 MeV, while tuk ~ 494 MeV. This 
explains why Figs.[16]and [T7]show remarkably small vari- 
ations of the kaon emitting source with increasing values 
of the transverse momenta of kaons. In case of protons, 
one would expect even smaller variations with increasing 
momentum, however, in this case the Monte Carlo sim- 
ulations indicate an increase of the effective source size 
with increasing momentum. Even in case of protons the 
dependence of the power-law exponent on the transverse 
momentum is negligible. Thus the transverse momentum 
dependence of the power-law exponent a is negligible in 
each considered case. 



D. Particle type dependence 

Particle type dependence of simulated source distribu- 
tions in various centrality classes and with a pair trans- 
verse momentum of 0.2 GeV/c < p t < 0.4 GeV/c is 
shown in Figs. [TIiriTil These source distributions all have 
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FIG. 20: Source distribution of pion, kaon and proton pairs 
with 0.2 GeV/c < p t < 0.4 GeV/c and 0-10% centrality. 
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FIG. 21: Source distribution of pion, kaon and proton pairs 
with 0.2 GeV/c < p t < 0.4 GeV/c and 30-50% centrality. 
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FIG. 22: Source distribution of pion, kaon and proton pairs 
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an approximate power-law tail, as expected in case of an 
anomalous diffusion effect. However, the power-law ex- 
ponents have rather different values for different particle 
types. This is well understandable as the total inelas- 
tic cross sections of these particles are rather different. 
Protons have the largest cross sections, pions the second 
largest, and kaons the smallest. Hence the protons have 
the shortest mean free path, the pions the second short- 
est, and the kaons have the largest mean free path at 
any given densities. This implies that the heaviest tail 
develops for kaons, the second heaviest tail for pions, and 
the proton source distribution is closest to the Gaussian 
distribution. However, it is rather difficult to make such 
a simple picture given that the HRC simulation includes 
a parameterized form of the strongly momentum depen- 
dent cross sections, when available, it uses Particle Data 
Group tables l37l. although not the most recent version 
of these tables [38J. Whenever PDG data were not avail- 
able, the HRC model utilized theoretical results on the 
total cross sections from Prakash, Prakash, Venugopalan 
and Welke [H. 



on heavy tails, i.e. it yields a power-law tail with an 
exponent, or Levy index of stability of a « 1.3 ± 0.1. De- 
tailed simulation results indicate that this exponent (as 
well as the tails of the S(r) relative coordinate distribu- 
tion) depends strongly on the particle cross sections or 
particle type, while for a given particle type it is remark- 
ably insensitive to centrality and momentum selections, 
which implies insensitivity to the actual number of colli- 
sions. The exponent value is rather different from a ~ 0.5 
predicted for the second order QCD phase transition in 
ref. [36[ . Also, at a second order QCD phase transition, 
the exponents are determined from universality class ar- 
guments, so in case of a second order phase transition, 
no particle type dependence is expected. Hence measur- 
ing the tails of particle production for protons and kaons 
seems to be a promising method to distinguish between 
a second order QCD phase transition and anomalous dif- 
fusion. 

The Levy distribution is called stable because the con- 
volution of two such distribution is also a Levy distribu- 
tion. In addition, the generalized central limit theorem 
says that the convolution of many (in the limiting case in- 
finitely many) elementary processes with the same (but 
arbitrary) probability distribution is a Levy stable dis- 
tribution. So it is natural, although still surprising that 
Levy stable, power-law tailed source distributions appear 
in the HRC model simulations. (The rescattering can be 
considered as such elementary process.) 

The extraction of the precise values of the power-law 
exponents, or the Levy index of stability, with reliable 
experimental statistical and systematic errors is an im- 
portant future task that allows experimental characteri- 
zation of anomalous diffusion or second order QCD phase 
transitions with a simple number. Measuring the trans- 
verse momentum, centrality, and particle type depen- 
dence of this exponent tells a lot about the production 
mechanism of particles in high energy heavy ion colli- 
sions. In particular, measuring the Levy index of stabil- 
ity for pions, kaons, and protons can serve as a promising 
experimental control possibility on the origin of the heavy 
tailed distribution. 



VII. CONCLUSIONS AND SUMMARY 

A power-law tail of the relative coordinate distribution 
S(r) appears in HRC model simulations for Au+Au col- 
lisions at RHIC energies. The simulation reproduces an 
important feature of recent PHENIX experimental data 
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